home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_6.2 / DPP / DPP.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  8KB  |  330 lines
  1. /* 
  2.  * dpp.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * D(raw)P(oint)P(attern)
  11.  *
  12.  * generate and draw random point, square or disk pattern
  13.  *
  14.  */
  15.  
  16. #include "dpp.h"
  17.  
  18. #define    ON            1
  19. #define    OFF            0
  20. #define    NEW_SEED    ON
  21. #define    SQR(a)        ( (a)*(a) )
  22.  
  23. /* globals */
  24. extern char *optarg;
  25. extern int optind, opterr;
  26. extern short tiffInput;         /* flag=0 if no ImageIn to set tags; else =1 */
  27.  
  28.  
  29. /*
  30.  * usage of routine
  31.  */
  32. void
  33. usage (char *progname)
  34. {
  35.   progname = last_bs (progname);
  36.   printf ("USAGE: %s outimg [-p t][-c f][-s x y] [-L]\n", progname);
  37.   printf ("\n%s generates and draws a random, ordered or clustered point pattern\n\n", progname);
  38.   printf ("ARGUMENTS:\n");
  39.   printf ("      outimg: output image filename (TIF format)\n\n");
  40.   printf ("OPTIONS:\n");
  41.   printf ("        -p t: select pattern type: R(andom), O(rdered), C(lustered)\n");
  42.   printf ("              t = R: random (default)\n");
  43.   printf ("              t = O: ordered\n");
  44.   printf ("              t = C: clustered\n");
  45.   printf ("        -c f: desired fractional coverage, 0<=f<=100\n");
  46.   printf ("      -s x y: output image size (default 512x512)\n");
  47.   printf ("          -L: print Software License for this module\n");
  48.   exit (1);
  49. }
  50.  
  51. /*
  52.  * function to draw points
  53.  */
  54. void
  55. draw_points (int *x, int *y, long n, struct Image *imgOut, int value)
  56. {
  57.   int ip;
  58.  
  59.   for (ip = 0; ip < n; ip++) {
  60.     setpixel (*(x + ip), *(y + ip), imgOut, value);
  61.   }
  62. }
  63.  
  64. /*
  65.  * function to generate a single random point coordinate
  66.  */
  67. Sp
  68. random_point (int nx, int ny)
  69. {
  70.   Sp pt;
  71.  
  72.   pt.x = (short) ((float) rand () * (float) nx / (float) RAND_MAX);
  73.   pt.y = (short) ((float) rand () * (float) ny / (float) RAND_MAX);
  74.   return (pt);
  75. }
  76.  
  77. /*
  78.  * function to generate random point coordinates
  79.  */
  80. void
  81. generate_random_points (int *xc, int *yc, long npix, int nx, int ny)
  82. {
  83.   int seed = 1;
  84.   int ip;
  85.   Sp pt;
  86.  
  87.   if (NEW_SEED == ON) {
  88.     printf ("...select new seed(enter 1 to reinitialize):");
  89.     scanf ("%d", &seed);
  90.   }
  91.  
  92.   srand (seed);                 /* must be called only once per series */
  93.   for (ip = 0; ip < npix; ip++) {
  94. #ifdef DEBUG
  95.     if (ip % 100 == 0)
  96.       printf ("    ip = %d\n", ip);
  97. #endif
  98.     pt = random_point (nx, ny);
  99.     *(xc + ip) = (int) pt.x;
  100.     *(yc + ip) = (int) pt.y;
  101.   }
  102. }
  103.  
  104. /*
  105.  * function to generate ordered point coordinates
  106.  */
  107. void
  108. generate_ordered_points (int *xc, int *yc, long npix, int nx, int ny, float coverage)
  109. {
  110.   int sepx;
  111.   int sepy;
  112.   int ix;
  113.   int iy;
  114.   int xnpts;
  115.   int ynpts;
  116.   int i;
  117.  
  118.   if (!npix)
  119.     return;
  120.   xnpts = (int) ((float) nx * (float) sqrt ((double) coverage));
  121.   ynpts = (int) ((float) ny * (float) sqrt ((double) coverage));
  122.   sepx = (nx / xnpts) + 1;
  123.   sepy = (ny / ynpts) + 1;
  124.   i = 0;
  125.   for (ix = 0; ix < nx; ix += sepx) {
  126.     for (iy = 0; iy < ny; iy += sepy) {
  127.       if (i >= npix)
  128.         break;
  129.       *(xc + i) = ix;
  130.       *(yc + i) = iy;
  131.       i++;
  132.     }
  133.   }
  134. }
  135.  
  136. /*
  137.  * function to generate clustered point coordinates
  138.  */
  139. void
  140. generate_clustered_points (int *xc, int *yc, long npix, int nx, int ny, long nxny)
  141. {
  142.   char buf[256];
  143.   long Nc;
  144.   int kc;
  145.   int ip;
  146.   int seed = 1;
  147.   Sp cur_seed_pt;
  148.   Sp new_pt;
  149.   int xsd, ysd;
  150.   int ppc;
  151.   double sigma;
  152.   double fxc, fyc;
  153.   double fac, rsq;
  154.  
  155.  
  156.   for (;;) {
  157.     printf ("How many clusters? ");
  158.     readlin (buf);
  159.     if (!(Nc = atol (buf)))
  160.       printf ("Sorry, invalid input, try again\n");
  161.     else
  162.       break;
  163.   }
  164.  
  165.   for (;;) {
  166.     printf ("Select SIGMA to set cluster width (width  proportional to NX/SIGMA): ");
  167.     readlin (buf);
  168.     if (!(sigma = (float) atof (buf)))
  169.       printf ("Sorry, invalid input, try again\n");
  170.     else
  171.       break;
  172.   }
  173.  
  174.  
  175.   if (NEW_SEED == ON) {
  176.     printf ("...select new seed(enter 1 to reinitialize):");
  177.     readlin (buf);
  178.     sscanf (buf, "%d", &seed);
  179.   }
  180.   srand (seed);
  181.  
  182.   ppc = npix / Nc;
  183.   ip = 0;
  184.   for (kc = 0; kc < Nc; kc++) {
  185.  
  186.  
  187.     cur_seed_pt = random_point (nx, ny);
  188.  
  189. /*
  190.  * ** restrict seed coordinate to interior of available AOI
  191.  * ** (to reduce chance of generating a cluster overlapping edge of AOI)
  192.  */
  193.     xsd = (int) ((nx / 2) + (2.0 * cur_seed_pt.x / nx - 1.0) * (nx / 8));
  194.     ysd = (int) ((ny / 2) + (2.0 * cur_seed_pt.y / ny - 1.0) * (ny / 8));
  195.     printf ("\ncluster: %d - seed point: %3d, %3d\n", kc, xsd, ysd);
  196.  
  197.     do {
  198.       new_pt = random_point (nx, ny);
  199.       fxc = (double) (2.0 * new_pt.x / nx - 1.0);
  200.       fyc = (double) (2.0 * new_pt.y / ny - 1.0);
  201.  
  202.       if ((rsq = SQR (fxc) + SQR (fyc)) < 1.0) {
  203.         fac = sqrt (-2.0 * log (rsq) / rsq);
  204.         *(xc + ip) = (int) (xsd + fac * fxc * (double) nx / (2 * sigma));
  205.         *(yc + ip) = (int) (ysd + fac * fyc * (double) ny / (2 * sigma));
  206.         //printf("\n%3d -  (%3d, %3d)\n", ip, *(xc+ip), *(yc+ip));
  207.         ip++;
  208.       }
  209.     } while (ip < ppc * (kc + 1));
  210.   }
  211.  
  212. }
  213.  
  214. void
  215. main (int argc, char **argv)
  216. {
  217.  
  218.   int i_arg;
  219.   char type = 'R';              /* type of pattern drawn */
  220.   int *xc, *yc;
  221.   float coverage;
  222.   long nxny;
  223.   long npix;
  224.   int nx = HEIGHT;
  225.   int ny = WIDTH;
  226.   struct Image *imgOut;
  227.  
  228. /* 
  229.  * cmd line options
  230.  */
  231.   static char *optstring = "p:c:s:L";
  232.  
  233.   nxny = nx * ny;
  234.   coverage = (float) (5.0);
  235.  
  236.  
  237. /*
  238.  * parse command line
  239.  */
  240.   optind = 2;
  241.   opterr = ON;                  /* give error messages */
  242.  
  243.   if (argc < 2)
  244.     usage (argv[0]);
  245.  
  246.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  247.     switch (i_arg) {
  248.     case 'c':
  249.       coverage = (float) atof (optarg);
  250.       break;
  251.     case 'p':
  252.       type = optarg[0];
  253.       break;
  254.     case 's':
  255.       if (!sscanf (argv[optind - 1], "%d", &nx) || !sscanf (argv[optind], "%d", &ny)) {
  256.         printf ("Error getting values for image size\n");
  257.         printf ("Will use defaults\n");
  258.         nx = WIDTH;
  259.         ny = HEIGHT;
  260.       }
  261.       else
  262.         optind += 1;
  263.       nxny = nx * ny;
  264.       break;
  265.     case 'L':
  266.       print_sos_lic ();
  267.       exit (0);
  268.     default:
  269.       printf ("...unknown command line option encountered\n");
  270.       exit (1);
  271.       break;
  272.     }
  273.   }
  274.  
  275.   npix = (long) (coverage * nxny / 100.0);
  276.   printf ("Image size for %s = %d x %d pixels\n", argv[1], nx, ny);
  277.   printf ("Coverage = %f percent of area, # pixels to be drawn = %ld\n", coverage, npix);
  278.   printf ("pattern = %c\n", type);
  279.  
  280. /*
  281.  * initialize graphics
  282.  */
  283.   imgOut = ImageAlloc ((long) ny, (long) nx, BPS);
  284.  
  285. /* reset tiffInput so that we write a grayscale file (i.e tags are not copied) */
  286.   tiffInput = 0;
  287.  
  288.   if ((xc = (int *) calloc ((size_t) npix, sizeof (int))) == NULL) {
  289.     printf ("\n...mem alloc for xc failed\n");
  290.     exit (1);
  291.   }
  292.   if ((yc = (int *) calloc ((size_t) npix, sizeof (int))) == NULL) {
  293.     printf ("\n...mem alloc for yc failed\n");
  294.     exit (1);
  295.   }
  296.  
  297.   switch (type) {
  298.   case 'R':                    /* Random point pattern */
  299.   case 'r':
  300.     generate_random_points (xc, yc, npix, nx, ny);
  301.     draw_points (xc, yc, npix, imgOut, WHITE);
  302.     break;
  303.   case 'O':                    /* Ordered point pattern */
  304.   case 'o':
  305.     generate_ordered_points (xc, yc, npix, nx, ny, coverage / (float) 100.0);
  306.     draw_points (xc, yc, npix, imgOut, WHITE);
  307.     break;
  308.   case 'C':                    /* Clustered point pattern */
  309.   case 'c':
  310.     generate_clustered_points (xc, yc, npix, nx, ny, nxny);
  311.     draw_points (xc, yc, npix, imgOut, WHITE);
  312.     break;
  313.   default:
  314.     printf ("...unknown pattern type: %c\n", type);
  315.     exit (1);
  316.     break;
  317.   }
  318.  
  319.  
  320. /*
  321.  * free memory
  322.  */
  323.   free (xc);
  324.   free (yc);
  325.  
  326.   printf ("\n...writing output file %s\n", argv[1]);
  327.   ImageOut (argv[1], imgOut);
  328.  
  329. }
  330.